home *** CD-ROM | disk | FTP | other *** search
/ Mac100% 1998 November / MAC100-1998-11.ISO.7z / MAC100-1998-11.ISO / ｵﾝﾗｲﾝｿﾌﾄ定点観測 / ﾕｰﾃｨﾘﾃｨ / Mops 3.2.sea / Mops 3.2 / Mops source / PPC source / test2 < prev    next >
Text File  |  1998-04-20  |  3KB  |  137 lines
  1. ¥ expint     Real Exponential Integral         ACM Algorithm #20
  2.  
  3. ¥ Forth Scientific Library Algorithm #1
  4.  
  5. ¥ Evaluates the Real Exponential Integral,
  6. ¥     E1(x) = - Ei(-x) =   int_x^¥infty exp^{-u}/u du      for x > 0
  7. ¥ using a rational approximation
  8.  
  9. ¥ This code conforms with ANS requiring:
  10. ¥      1. The Floating-Point word set
  11. ¥      2. The immediate word '%' which takes the next token
  12. ¥         and converts it to a floating-point literal
  13. ¥
  14.  
  15. ¥ Collected Algorithms from ACM, Volume 1 Algorithms 1-220,
  16. ¥ 1980; Association for Computing Machinery Inc., New York,
  17. ¥ ISBN 0-89791-017-6
  18.  
  19. ¥ (c) Copyright 1994 Everett F. Carter.  Permission is granted by the
  20. ¥ author to use this software for any application provided the
  21. ¥ copyright notice is preserved.
  22.  
  23.  
  24. CR .( EXPINT     V1.1                  21 September 1994   EFC )
  25.  
  26. variable v1
  27. variable v2
  28. variable v3
  29. variable v4
  30. variable v5
  31. variable v6
  32. variable v7
  33. variable v11
  34. variable v12
  35. variable v13
  36. variable v14
  37. variable v15
  38. variable v16
  39. variable v17
  40. variable v18
  41. variable v21
  42. variable v22
  43. variable v23
  44. variable v24
  45. variable v25
  46. variable v26
  47. variable v27
  48. variable v28
  49.  
  50. : expint ( --, f: x -- expint[x] )
  51.         FDUP
  52. (*        % 1.0 F< IF
  53.                     FDUP % 0.00107857 F* % 0.00976004 F-
  54.                     FOVER F*
  55.                     % 0.05519968 F+
  56.                     FOVER F*
  57.                     % 0.24991055 F-
  58.                     FOVER F*
  59.                     % 0.99999193 F+
  60.                     FOVER F*
  61.                     % 0.57721566 F-
  62.                     FSWAP FLN F-
  63.                 ELSE
  64. *)
  65.                 FDUP v4 f@ F+
  66.                     FOVER F*
  67.                     v5 f@ F+
  68.                     FOVER F*
  69.                     v6 f@ F+
  70.        db
  71.                     FOVER F*
  72.                     v7 f@ F+
  73. (*        v11 f@
  74.         v12 f@
  75.         v13 f@
  76.         v14 f@
  77.         v15 f@
  78.         v16 f@
  79.         v17 f@
  80.         v18 f@
  81.         v21 f@
  82.         v22 f@
  83.         v23 f@
  84.         v24 f@
  85.         v25 f@
  86.         v26 f@
  87.         v27 f@
  88.         v28 f@
  89. *)
  90.                     FOVER
  91.                     FDUP v1 f@ F+
  92.                     FOVER F*
  93.                     v2 f@ F+
  94.                     FOVER F*
  95.                     v3 f@ F+
  96.                     FOVER F*
  97.                     % 3.9584969228 F+
  98.  
  99.                     FSWAP FDROP
  100.                     F/
  101.                     FOVER F/
  102.                     FSWAP % -1.0 F* FEXP
  103.                     F*
  104.  
  105.             ¥    THEN
  106. ;
  107.  
  108.  
  109. ¥ test code generates a small table of selected E1 values.
  110. ¥ most comparison values are from Abramowitz & Stegun,
  111. ¥ Handbook of Mathematical Functions, Table 5.1
  112.  
  113. : expint_test ( -- )
  114.  
  115.         CR
  116.         ."   x    E1(x) exact      ExpInt[x] " CR
  117.  
  118.  
  119.       ."  0.5   0.5597736      "
  120.       % 0.5 expint  F. CR
  121.  
  122.       ."  1.0   0.2193839      "
  123.       % 1.0 expint   F. CR
  124.  
  125.       ."  2.0   0.0489005      "
  126.       % 2.0 expint    F. CR
  127.  
  128.       ."  5.0   0.001148296    "
  129.       % 5.0 expint    F. CR
  130.  
  131.       ." 10.0   0.4156969e-5   "
  132.       % 10.0 expint   F. CR
  133.  
  134. ;
  135.  
  136.  
  137.